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ABSTRACT 

A  method  is  proposed  for  the  estimation  of  the  transfer 
function  of  a  linear,  time-invariant  system  with  no  numera- 
tor dynamics,  from  random  input  and  output  data.   The  method 
employed  utilizes  the  Fast  Fourier  Transform  Algorithm  and 
Least  Squares  Estimation  to  obtain  the  coefficients  of  the 
system's  transfer  function. 

The  procedure  has  been  modeled  in  FORTRAN  IV  on  an 
IBM-360  computer.   The  results  of  simulation  show  the 
feasibility  of  estimating  the  order  of  the  transfer  function 
and  its  coefficients. 
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I.   INTRODUCTION 

The  problem  of  identifying  a  system  (or  plant)  has  been 
the  object  of  much  thought  and  effort  as  evidenced  by  the 
Identification  in  Automatic  Control  Systems  Symposium  [1]. 
The  effort  has  been  primarily  expended  on  methods  which  are 
suitable  to  a  time  domain  analysis.   Identification  in  the 
frequency  domain  has  in  general  involved  lengthy  computa- 
tions to  obtain  a  Fourier  Transform  if  the  computations  are 
done  digitally,  or  a  tradeoff  between  many  parallel  narrow 
band  filters  or  repetitive  playback  of  the  system  input 
and  output  if  the  analysis  is  done  by  analog  methods.   The 
method  to  be  used  in  this  investigation  has  been  constrained 
to  make  use  of  the  Fast  Fourier  Transform  algorithm.   There- 
fore digital  computations  and  techniques  logically  follow. 
The  Fast  Fourier  Transform  (FFT)  algorithm  reduces  the 
number  of  complex  multiplications  required  in  the  evaluation 
of  the  Fourier  Transform  and  thus  may  provide  a  means  of 
performing  system  identification  on-line  or  approximately  in 
real  time.   The  FFT  is  discussed  in  Appendix  A. 


II.   PROBLEM  STATEMENT 

The  mathematical  model  of  a  linear,  time-invariant 
system  may  be  expressed  as  a  differential  equation  or  as 
a  transfer  function.   The  determination  of  the  model 
parameters  from  observations  on  the  system  which  is  to  be 
modeled  is  generally  called  an  identification  problem. 
This  investigation  considers  the  modeling  of  a  system, 
which  is  characterized  by  a  vector  differential  equation, 
such  that  the  coefficients  of  the  equivalent  transfer 
function  can  be  estimated. 

Let  a  linear  system  transfer  function  be  written  as 


wry  ■  ditt  >    D(f)  -    *    d-(J")n  • 


n=0 


n 


(1) 


rth 


The  equivalent  N    order  differential  equation  is  written 
N 


,  d  x( t )     / ,  n 

L      d — -  =  v(t)  . 

n   n   ,,n 
n=0      dt 


(2) 


dN  =  1.0 

From  (2)  the  equivalent  N-vector  differential  equation 
(3)  can  be  obtained. 
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x(t)  =  Ax(t)  +  Bv(t) 
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N-l    N 


x,T(t)    1 


v(t) 
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The  solution  of  (3)  is 


t 


x(t)  =  <D(t,t0)x(t0)  +  /   *(t,T)  •  B  •  v(T)di,     (4) 

t0 

where  x(t)  is  the  system  state  vector  and  the  element  x.,(t) 
is  the  system  output  x(t).   The  discrete  version  of  (4), 
which  is  derived  in  Appendix  B,  is  used  as  the  equivalent 

representation  of  (1)  and  the  coefficients  ctL,  d,  ,  d? d^ 

are  estimated.   The  estimation  scheme  is  represented  by 
Figure  1. 

The  investigation  is  to  be  constrained  so  that  v(t)  is 
broad-band  white  noise  and  the  Fast  Fourier  Transform 
algorithm  must  be  used  in  the  identification  scheme. 
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UNKNOWN  SYSTEM  (EQUATION  4) 
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Figure  1.   Estimation  Scheme  Representation 
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III.   SYSTEM  MODEL  FOR  A  DIGITAL  COMPUTER 

A  linear,  time-invariant  system  represented  by  the 
vector  differential  equation  (5)  has  equation  (6)  as  its 
solution. 

x(t)  =  Ax(t)  +  Bv(t)  (5) 


A(t-t  )  t   fl/+.  t) 

X(t)   =  E        U  X(tn)   +  /    £Mt"T;Bv(T)dT  (6) 

—   u     t 


If  v(t)  is  deterministic  and  the  convolution  integral 

can  be  solved,  then  a  discrete  solution  of  (5)  can  be 

A(t-tQ) 
modeled  on  a  digital  computer.  e  can  be  converted 

to  a  square  matrix  by  using  the  Laplace  transform  method 

[2,  pg.  316].   However,  since  v(t)  is  not  deterministic  the 

convolution  integral  cannot  be  solved,  so  an  approximation 

to  (6)  is  required. 

Considering  only  the  unforced  response  of  (6), 

A(t-tn) 
x(t)  =  e  U  x(tQ)  .  (7) 

At  some  instant  of  time  t  =  tM  =  t„  +  MT, 

x(tM)  =  £AMT  x(tQ)  (8) 


which  is  identical  to  the  unforced  response  of  equation  (3) 
of  Appendix  B  where  x(MT)  =  x(tQ  +  MT)  =  x(tM), 
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AT 
x((m+l)T)  =  e   x(mT)  (9) 

x(T')  =  eATx(0)  =  eATx(t0) 

x(2T)  =  eATx(T)  =  eA2Tx(t0) 


x(MT)  =  eAMTx(tQ)  (10) 

Therefore,  since  (10)  is  equivalent  to  (8)  which  is  the 

discrete  solution  of  (7),  then  (9)  is  the  discrete  solution 

for  the  unforced  response  of  (6)  for  t  =  t,     M\=  tn+  (m+l)T 

^  (m+1)    0 

or  in  the  notation  of  Appendix  B,  t  =  (m+l)T. 
The  convolution  integral  from  (6), 


;(t)  =  /  £A(t-T)Bv(T)dT,  (11) 

t0 


cannot  be  forced  to  yield  an  exact  discrete  solution  without 
performing  the  indicated  integration.   Therefore,  the  sam- 
pled data  approach  in  Appendix  B  is  necessary.   Normally 
this  approach  yields 


g((m+l)T)  =  /   eA(T  T)Bv(mT)dT,  (12) 


where  v(mT)  is  assumed  to  be  constant  for  mT^t<(m+l)T 
[2,  pg.  3^0].   This  approach  presupposes  that  the  bandwidth 
of  v(t)  is  small  compared  to  ^.   Since  v(t)  is  assumed  to 
be  broadband  noise,  T  should  be  small.   However,  as  T 
decreases,  the  number  of  iterations  required  in  the  simula- 
tion for  a  given  time  span  of  x(t)  increases.   Each  iteration 

14 


requires  that  both  the  unforced  and  forced  responses  of 
equation  (6)  be  computed.   T  should  then  be  large  to  min- 
imize computation  and  small  to  allow  the  maximum  transfer 
of  information  to  the  system  concerning  the  variations  of 
v(t).   Equation  (12)  of  Appendix  B  was  chosen  as  a  compro- 
mise.  By  using  this  equation,  v(t)  can  be  sampled  at  a 
rate  greater  than  ~   and  the  unforced  response  of  equation 
(6)  need  only  be  computed  once  every  T  time  interval. 


15 


IV.   SYSTEM  IDENTIFICATION  SCHEME 

A  continuous,  time-invariant  system  is  to  be  identified 
by  estimating  the  coefficients  of  its  transfer  function. 
The  input  to  the  system  is  broad-band  white  noise  of  known 
mean  and  variance.   The  identification  scheme  is  depicted 
by  Figure  2,  where  F  denotes  the  Fourier  Transform,  i.e., 

V(f)  =  F[v(t)] 

X(f)  =  F[x(t)] 

v(t)  -  broadband  white  noise 

x(t)  -  system  output 

The  object  of  this  investigation  is  to  make  use  of 
the  Fast  Fourier  Transform  (FFT,  Appendix  A)  in  place  of  F, 
and  Least  Squares  Estimation  in  the  computations  for  the 
transfer  function  coefficients. 


NOISE 
GENERATOR 


vCt] 


I 


UNKNOWN 
SYSTEM 


zlll 


COMPUTATIONS  FOR 
TRANSFER  FUNCTION 
IDENTIFICATION 


I 


TRANSFER  FUNCTION  COEFFICIENTS 


Figure  2.   Identification  Block  Diagram 
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In  Figure  2,  v(t)  is  the  output  of  a  random  process. 
In  particular  it  is  the  output  of  a  gaussian  random  process 
Since  v(t)  is  random,  then  x(t),  V(f)  and  X(f)  are  random 
processes.   V(f)  and  X(f)  are  not  themselves  deterministic; 
however,  from  Goldman,  [3,  pg.  279] 

H(f)  =  X(f)/V(f)  ,  (13) 

where  H(f)  is  the  frequency  response  of  the  unknown  system. 
Since  H(f)  is  assumed  to  be  of  the  form 


H(f)  = 


N  D(f ) 

E   d(j27Tf)n 
n=0   n 

dN  =  1.0, 
then  (13)  can  be  written  as  [Appendix  C,  equation  (2)] 

X(f)  •  D(f)  =  V(f) .  (14) 

From  (14)  a  residue  r(f)  is  formed, 

r(f)  =  X(f)  •  D(f)  -  V(f).  (15) 

Through  the  Fourier  Transform,  measurements  on  x(t)  and 
v(t)  yield  X(f)  and  V(f).   A  modified  least  squares  esti- 
mation procedure  then  uses  r(f)  from  (15)  and  its  complex 
conjugate  to  form 

F 
Q  =  /   r(f)  •  r*(f)df  (16) 

0 
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where  r*(f)  is  the  complex  conjugate  of  r(f).   The  object 
is  then  to  minimize  Q  with  respect  to  the  d  's,  the 
coefficients  of  D(f ) . 

Since  the  FPT  [Appendix  A]  is  used  in  place  of  the 
Fourier  Transform  of  Figure  2  and  the  unknown  system  is 
modeled  in  accordance  with  Appendix  B,  the  identification 
scheme  takes  on  the  form  of  Figure  3- 

Since  u(t)  is  broadband  noise,  the  low-pass  filter  is 
a  necessary  addition  to  the  identification  scheme  to  prevent 
foldover  of  the  alias  frequencies  of  V. .   The  low-pass 
filter  used  is  a  discrete  version  of  an  R-C  filter  with 
its  break  point  on  a  Bode  plot  being  greater  than  the 
expected  bandwidth  of  the  system  under  investigation.   The 
sampling  rate  of  the  inputs  to  the  FFT ' s  was  then  chosen 
to  be  at  least  twice  as  great  as  the  low-pass  filter  cut- 
off frequency,  that  frequency  at  which  the  filter  output 
drops  by  3db . 

X.  and  V.  are  the  results  of  the  Discrete  Fourier  Trans- 
1       l 

forms  of  x(t)  and  u(t)  which  are  observed  through  a  time 
window  of  t  seconds,  i.e. 

X.  =  X(f . )  =  X(-) 

1         1         T 

V.  =  V(f . )  =  V(-)  , 
1        1        T 

where  t  is  the  time  interval  over  which  u(t)  and  x(t)  are 
observed.   As  the  index  i  goes  from  -tW  to  tW,  X.  and  V. 
represent  the  complex  values  of  the  Fourier  Transforms  of 
x(t)  and  u(t)  at  the  discrete  frequencies  -W  to  +W. 
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The  X.'s  and  V.'s  exist  only  at  frequencies  which  are 

multiples  of  — .   If  x(t)  and  v(t)  contain  only  frequencies 

which  are  exact  multiples  of  —  then  the  X.'s  and  V.'s 

^        T  11 

will  be  the  exact  Fourier  coefficients  of  x(t)  and  v(t)[6]. 
However,  since  it  Is  unlikely  that  the  only  frequencies 
present  are  multiples  of  — ,  these  other  frequency  components 

will  cause  the  X.'s  and  V.'s  to  be  altered.   It  would  be 

1        1 

advisable  to  weight  the  frequency  spectra  of  x(t)  and 
v(t)  with  the  hamming  or  hanning  functions  [4,  5,  6,  Appen- 
dix A]  to  minimize  the  effects  of  these  other  frequency 
components.   This  weighting  was  not  performed,  although  it 
is  admittedly  advisable,  since  computer  computation  time 
would  be  increased.   The  trade-off  of  less  computation  time 
for  increased  accuracy  was  felt  to  be  justified  since  much 
of  the  computer  time  used  in  this  investigation  had  to  be 
charged  to  other  projects. 
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V.   RESULTS 

Many  of  the  variables  in  the  computer  simulation  of  the 
"unknown"  plant  and  the  identification  scheme  [Table  Dl, 
Appendix  D]  are  peculiar  to  the  approach  taken  in  this 
investigation.   For  an  actual  application  of  the  identifi- 
cation scheme  the  "unknown"  plant  is  not  simulated  and  only 
the  sampling  rate  and  number  of  samples  to  be  taken  must 
be  determined.   The  expected  or  known  signal  bandwidth 
dictates  the  sampling  rate,  where  the  minimum  rate  is  the 
Nyquist  sampling  rate.   An  increase  in  the  number  of 
samples  taken  increases  the  frequency  resolution  of  the 
FFT  output  and  the  accuracy  of  the  estimated  transfer 
function  coefficients.   However,  it  also  increases  compu- 
tation time  and  computer  memory  required. 

A  tenth  order  system  such  as  equation  (17)  was  assumed 
for  each  trial. 


X(S)       1.0  M7s 

vTsT  "  To (17) 

n 

Z   d  S 

n    n 

n=0 


The  low-pass  filter  shown  in  Figure  3  is  represented  by 
equation  (18) . 

V(S)    _      2tt    •    100  nfn 

UTST        S  +  2tt  •  100  {       } 

Table  1  presents  the  most  encouraging  of  the  results 
of  this  investigation.   Trials  1  and  2  represent  first-order 
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systems  with  real  S-plane  poles.   Trials  3  and  4  represent 
second-order  systems  with  real  and  complex  poles  respect- 
ively.  Trial  5  represents  a  second-order  system  with 
complex  poles  and  trial  6  is  a  third-order  system  with 
real  poles  of  -10,  -20  and  -30.   Only  the  five  lowest-order 
coefficients  of  the  estimates  are  given  in  Table  1  since 
the  higher-order  estimates  become  negligibly  small.   The 
results  given  are  the  averages  of  fifty  consecutive  esti- 
mates made  for  each  system's  coefficients. 
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TABLE  1 


ESTIMATION  RESULTS 
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VI.   CONCLUSIONS 

The  identification  scheme  presented  in  this  investiga- 
tion does  provide  good  estimates  of  the  transfer  function 
coefficients  for  an  unknown  system  with  no  numerator 
dynamics.   The  estimated  coefficients  in  Table  1  are  seen 
to  be  within  50%  of  the  actual  coefficients  for  all  cases 
except  when  the  actual  coefficient  equals  zero.   For  these 
coefficients  the  estimates  rapidly  approach  zero  since  the 
estimates  are  for  higher  order  coefficients  than  those  that 
exist  In  the  actual  system.   The  identification  scheme 
therefore  provides  an  indication  of  the  order  of  the  system. 

The  need  for  a  low  pass  filter  on  the  input  to  the 
"unknown"  system  has  not  been  substantiated.   Trials  1 
through  6  of  Table  1  indicate  that  the  filter  can  be  deleted 
with  no  disastrous  effects. 

The  feasibility  of  using  the  Past  Fourier  Transform  as 
a  tool  for  system  (plant)  identification  has  been  shown. 
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VII.   RECOMMENDATIONS 

This  investigation  has  been  unduly  long  and  involved 
due  to  the  requirement  that  all  the  modeling  be  done  on  a 
digital  computer.   Many  of  the  parameters  listed  in  Table 
Dl  of  Appendix  D  required  a  large  number  of  trial  runs  to 
determine  a  best  set  of  values  for  the  overall  simulation. 
These  trials  were  necessary  to  aid  in  finding  the  minimum 
iteration  rate  and  input  signal  sampling  rate  required  to 
adequately  simulate  the  continuous  "unknown"  plant  since 
these  rates  are  fundamental  to  the  total  simulation  time 
(not  to  be  confused  with  identification  time).   It  is 
suggested  that  any  further  investigation  in  this  area  be 
done  on  a  Hybrid  computer  so  that  the  continuous  plant  may 
be  run  on  the  analog  portion  of  the  computer  and  the  iden- 
tification scheme  on  the  digital  portion. 

As  a  starting  point  the  identification  scheme  requires 
that  a  guess  be  made  as  to  the  maximum  number  of  coeffi- 
cients in  the  unknown  plant's  transfer  function.   A  major 
effect  of  this  guess  is  that  it  determines  the  dimensions 
of  matrix  W  of  equation  (15),  Appendix  C.   Since  the  inverse 
of  W  is  required  it  would  appear  that  its  dimensions  should 
be  kept  as  small  as  possible.   However,  it  has  been  observed 
that  the  dimensions  of  W  that  produce  the  most  nearly  cor- 
rect estimates  of  the  transfer  function  coefficients  are 
not  the  smallest  (consistent  with  the  number  of  coefficients 
known  to  exist)  nor  the  largest  (which  depend  on  the  amount 
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of  computer  memory  reserved  for  W) .   The  best  dimensions 
for  W  are  different  even  for  different  plants  of  the  same 
order.   It  is  suggested  that  any  future  investigation  util- 
izing this  identification  scheme  should  involve  a  method 
of  avoiding  the  inverse  of  W  or  attempt  an  analytic  solution 
of  W"1. 
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APPENDIX  A 
The  Fast  Fourier  Transform  (FFT) 

The  FFT  is  used  to  evaluate  the  complex  Fourier  trans- 
form of  a  sequence  of  complex  numbers.   It  is  a  finitie 
discrete  Fourier  transform  (DFT)  with  the  restriction  that 
the  number  of  samples  in  the  complex  sequence  is  constrained 
to  be  an  exact  power  of  two.   The  power-of-two  constraint 
on  the  number  of  samples  allows  implementation  of  an  algor- 
ithm on  a  digital  computer  for  rapid  computation  of  the 
finite  DFT. 

The  mechanics  of  the  FFT  and  Its  characteristics  has 
been  widely  discussed  in  the  literature  [7,  8,  9,  10,  11, 
12,  13].   A  few  of  the  characteristics  will  be  presented 
here  . 

From  Goldman  [3] 3    a  function  g(t)  which  exists  only  for 

— T      T 

^p-  < t  <  —  and  which  has  a  Fourier  transform  G(f)  which  exists 

only  for  -W  < f  < W  can  be  expressed  by  the  DFT  pair 

/  m  v    1   ™    n  /  n  v   J°  ~2TW~  f  n  v 

n=-TW 

TW  4  2  mn 

n,ns         1   v     /  m  x   J  2TW  (0s 

G(T}  =  2W   *    g(2W}  £  (2) 

m=-TW 


where    S(^)  =  g(t),  for  t  =  ^g 


G(£)  =  G(f),  for  f  -  ^ 


28 


g(t)  and  G(f)  are  completely  determined  by  their  2TW+1 
sample  values  g(py)  and  G(~)    respectively. 

The  assumptions  leading  to  (1)  and  (2),  that  g(t)  exist 

-T      T 
only  for  -x-  <  t  <  p-  and  G(  f )  exists  only  for  -W  <  f  <  W,  are 

physically  unrealizable.   For  example,  consider  the  case 

where  g(t)  is  a  section  of  a  continuous  function  x(t)  such 

that 

x(t)  is  continuous  for  -°°  <  t  <  °°  (3) 

g(t)  =  x(t)  •  y(t) 

y(t)  =  u(t+|)  -  u(t-|) 

u(t)  is  the  unit  step  function  . 

oo 

Then     G(f)  =  /   X(q)Y(f-q)dq  (4) 

_oo 

and     Y(f)  =  T  sll\<;fT)  .  (5) 

Y(f)  and  therefore  G(f)  are  continuous  for  -°°  <  f  <  °°  . 

By  observing  (sampling)  a  section  of  x(t)  such  that 

-T     T 
g(t)  =  x(t)  for  -p- <  t  <—  the  assumption  leading  to  (1) 

and  (2)  are  violated. 

If  a  y(t)  in  (3)  is  chosen  so  that  g(t)  =  x(t)-  y(t) 

— T      T 
has  negligible  value  outside  of  the  range  -x-  <  t  <  x-   and  G(f) 

is  negligible  outside  the  range  -W  < f <  W,  then  (1)  and  (2) 

can  be  made  to  be  approximately  exact.   y(t)'s  to  satisfy 

the  above  include  the  Hanning  and  Hamming  functions  [4,  51 

and  the  Parzen  spectral  window  [14]. 
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For  the  FFT,  (1)  and  (2)  are  normally  expressed  as 

N-l       J  ~^~ 
g(m)  =   Z   G(n)e  (6) 

n=0 

l  N_1  ~J  ~w~ 

G(n)  =  ±  E   g(m)  e        ,  (7) 

m=0 


where    g(m)  =  gC^y) 

G(n)  =  G(£) 

n,m  =  0  to  N-l 
N  =  2k  . 

Some  authors  put  the  weighting  jt  on  (6)  rather  than  (7) 
However,  the  FFT  used  in  this  investigation  is  from  the 
IBM/360  Scientific  Subroutine  package  ( 360A-CM-03X)  and  in 
this  algorithm  the  rr   weighting  is  placed  on  (7). 
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APPENDIX  B 

Discrete  Solution  for  x(t) 

Given  the  vector  differential  equation  (1)  it  is  desired 
to  obtain  a  discrete  solution  for  x(t). 

x(t)  =  Ax(t)  +  Bv(t) ,  (1) 

T 
where    x(t)  =  [x-,(t)  Xp(t)  xN(t)] 

T 
x(t)  =  [x1(t)  x2(t)  ^N(t)] 

A  is  an  (N,N)  constant  matrix 
B  is  an  (N,l)  constant  vector 
u(t)  is  a  scalar  . 

The  superscript  T  denotes  the  transpose  of  the  vector.   Prom 
Ogata  [2]  the  continuous  time  solution  of  (1)  is 

A(  t-tn)  t  i\(f-_T\ 

x(t)  =  e      U   x(t  )  +  /  eR{Z    T)    Bv(t)  di,   (2) 

f 

A(t-tQ) 
where  e        is  the  transition  matrix  $(t,tn) 

For  a  discrete  solution  of  (2),  where  x(t)  is  sampled 
at  intervals  of  time  T,  let 

tQ  =  nT 
t  =  (n  + 1)T  . 
Then 

x  [(n+l)T]  =  £ATx(nT)  +  /        eA  [  (n+1 )T"T ]Bv( t )dx   (3) 

nT 
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The  usual  procedure  is  to  assume  that  v(t)  =  v(nT)  for, 
nT  <.  t  <(n+l)T.   The  implementation  of  v(t)  in  this  inves- 
tigation possesses  characteristics  which  make  it  advisable 
to  sample  this  forcing  function  at  as  high  a  rate  as 
possible.   Since  v(t)  is  assumed  to  be  broadband  noise,  it 
is  varying  very  rapidly  with  respect  to  time.   Therefore  T 
should  be  very  small  if  v(t)  is  to  be  fairly  constant  over 
the  sampling  interval.   Also  it  must  be  noted  that  in  the 
simulation  of  the  identification  problem  v(t)  is  represented 
by  a  sequence  of  numbers  from  a  random  number  generator. 
The  magnitude  of  the  numbers  in  this  sequence  are  indepen- 
dent of  T  in  that  if  five  numbers  are  produced  they  will 
have  the  same  magnitudes  whether  they  are  produced  at 
intervals  of  T  or  T/2 .   The  implication  then  is  that  for 
any  two  consecutive  numbers  from  the  generator,  such  as  v(i) 
and  v(i+l),  those  which  are  produced  at  T/2  second  intervals 
represent  a  v(t)  that  has  a  higher  rate  of  change  with 
respect  to  time  (and  a  higher  frequency)  than  those  that 
are  produced  at  T  second  intervals,  i.e., 


v(i+l)-v(i)  >  v(i+l)-v(i) 


\  T/2  T       J 

Therefore  if  the  v's  in  (3)  can  be  generated  at  intervals 
smaller  than  T,  then  the  effective  bandwidth  of  the  noise 
can  be  made  wider  as  the  sampling  interval  of  v(t)  becomes 
smaller.   Some  averaging  of  the  rapid  fluctuations  of  v(t), 
weighted  by  the  system  parameters  in  the  matrix  A,  will  be 


32 


made  by  the  integral  in  (3).   The  result  is  that  v(t)  can 
be  sampled  at  a  high  rate  (higher  than  „)  and  the  unforced 
response  in  (3)  need  only  be  computed  every  T  seconds.   To 
implement  this  approach  let, 


(m+1)  )ry 

M-l      M  LlJ± 

x[(n+l)T]  =  eATx(nT)+  I    /  eA [ (n+1 )T-t ] Bv( x )dx 

m=°  (^+n)T  (4) 


v(t)  =  v[(g+n)T],  (g+n)T  ±T<(^+n)T  . 

Since  B  is  a  constant  vector  and  v(t)  will  be  a  constant 
over  the  interval  of  integration, 


(m+1in)T 
x[(n+l)T]  =  £Aix(nT)+  Z         f  eA[(n+l  )T-t  JdT-B-v  [(— i-n)T]  . 

m=0  fULfn^T 

(M+n)T  (5) 

To  remove  the  parameter  n  from  the  integral  in  (5)  let 


t    =    nT-A 
A    =    nT-i 

Then 

_(Hi+l)T 

M-l  M 

x[(n+l)T]=£;ATx(nT)  +  £AT      Z    -/  eAAdA  •  B  •  v  [  (jj|+n )T]  . 

m=0   -mT 

M  (6) 

To  remove  the  negative  limits  on  the  integral  in  (6)  let 

mT 
M 

mT 
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X  =  "T-  M 


t  -  -A- 


Then  (6)  becomes 

T 

—         mT 

M-1   M   -A(t+— ) 
x[(n+l)T]  =  £Mx(nT)  +  EM   £   /   e      '   dx -B- v[ (g+n)T] .  (7) 

m=0  0 


Since  the  integral  in  (7)  is  not  a  function  of  m,  then 


M        M-1 
x[(n+l)T]  =  £ATx(nT)  +  £AT  /   e~ATdx   E   e~A  5^Bv[(g+N)T]  .   (8) 

0         m=0 


To  evaluate  the  integral  in  (8)  on  a  digital  computer, 

-At 
e     is  expanded  in  an  infinite  series. 

,    -At,     r         „   (-At)   , 
/   e    dx  =  /    1      - — .  ;    dT 

0  0   1=0    1- 

T 

»   ,  .a     M   . 
v   (-A)      r  1, 

I p —  •  /    T  dT 

i=0    1*     0 


=   Z 


T   (i+1) 
(-A)1      M 


i-o  1:      '  (1+1) 


,-AT    i 

oo             (  ) 

T      r       v    M    ; 

M    .;Q    (i+1)! 

T 

M 

(9) 


-At 
Equation  (9)  is  the  solution  of  /   e    dT  to  be  used.   How- 

0 
ever,  if  A  is  known  to  be  nonsingular,  then 

T 

M 

r         -AT,     .-1  rT   -AT,  Mnv 

J      e        dT  =  A    [I-e-ss— J  (10) 

0 
where  I  is  the  identity  matrix  of  the  same  dimensions  as  A. 
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A~   should  exist  except  for  those  cases  where  the  system 
characterized  by  (1)  has  a  pole  at  S=0.   However,  to  avoid 
compounding  errors  by  taking  an  unnecessary  matrix  inverse 
on  the  computer,  (9)  is  the  solution  of  the  integral  to  be 
used . 

From  (8)  and  (9)  and  letting  ^  =  ™, 


x[(n+l)T]  = 


£ATx(nT)+eAT-A-   Z 

i  =  0 


OR 


/  ., vi   M-l   .  , 

|^_.   E  £-AmA-B.v(nT+mX) 

(1+1)!   m=0 

(ID 


x[(n+l)T]  =  $(T)x(nT)  +  T (T,M) v(nT ,M) 


(12) 


where 


x(nT)  is  the  (N,l)  state  vector  of  the  plant  at 
time  t=nT 

AT 
$(T)  is  the  (N,N)  transition  matrix  e 

T  is  the  transition  time  between  the  system's  states 


v(nT,M)  is  the  (M,l)  control  vector 


v(nT) 
v(nT+X) 


v(nT+(M-l)X) 


r(.T,M)    is    an    (N,M)    transfer  matrix   such   that 


r(T,M)  =  $(T) 


.     _       (-AA)1 

\l0  mm 


• 

'    -AX*    -2AX* 

I  .  £             .  £                

-(M-l)AX     >B 
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APPENDIX  C 


Implementation  of  the  System  Identification  Solution 


A  given  system  is  represented  by  Figure  CI,  where 


K 


k 


D(oO  =  T,      d,  (joo)  .   It  is  desired  that  the  coefficients 


k=0 


k 


of  D(co)  be  identified  through  measurements  on  v(t)  and  x(t), 


v(t) 


V((D.) 


1 

x(t) 

U(w) 

X(u>) 

Figure  CI.   System  Representation  in  the 
Frequency  Domain 


Since 


1 
DT^T 


X(o>) 

vTwT  » 


(i) 


then 


X(oo)D(to)  =  V(co) 


(2) 


If  X(to)  and  V(w)  are  known  for  discrete  values  of  w  where 


oj  ->-  co.  then 
i 


X(o).  )D(co.  )  =  V(o).  )  . 
l     i        l 


(3) 


XC.to.)  andVCco.)  result  from  application  of  the  Fast 
Fourier  Transform  algorithm  [Appendix  A]  to  measurements 
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of  x(t)  and  v(t).   Therefore  X(w.)  and  V(w.)  are  assumed 
to  be  known.   If  D(co.)  in  (3)  is  replaced  by  D(w.  )  ,  the 
estimate  of  D(w.),  then  (3)  becomes 

X(o)i)D(o)1)  =  V(w1)  (4) 

Let  r(w.)  equal  the  residue,  or  error,  due  to  the 

estimate  D(co.  )  . 

i 

r(ooi)  =  X(u)1)D(o)1)  -  V(ok)  (5) 

Then  applying  a  least  squares  estimation  procedure  to 


I 
Q  =   £   r(oo.  )  •  r*(a), ),  (6) 

i  =  0     x        x 


*  denotes  the  complex  conjugate 
Q  is  to  be  minimized. 
From  (5)  and  (6 ) , 


[X(io.  )D(uk)-V(co.)]  •  [X(a3.  )D(oo.  )-V(oj.  )]* 


i  =  0 


(7) 


Let 


X((D.  )  =  X. 

l      l 


V(u>.  )  =  V. 
l      l 


D(co,  )  =   Z   d,  (jw,  )k  =  D. 
1     k=0   K    x       x 
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Then,  setting  the  partial  of  Q  with  respect  to  a  coeffi- 
cient  d   equal  to  zero, 

SO      Z      -  sS*i       d^i 

2$-   -   S   (X.D.-V.  )X*  -7T+  +    X,  -si  (X.D.-V,  )*  =  0,     (8) 

3d     1  =  0        x    x      x         ±    3(5  dd  x 

m  mm 


where 

3D, 

-^=  (j^)m  (9) 


3d 
m 


9D*.  Q    K  „        . 

-S —  =  -r-  ^  d,  (-jo),  )   =  (-1)  ( j  oo  )    .       (10) 

9d  3d   k=0  K                  1 

m  m 


From  (8),  (9)  and  (10), 


I      X,X*   [D,  (-l)m(Ju>,  )m  +  D*  (Jed.  )m] 

1  =  0 


I 

Z   V.X*.  (-l)m(jaj.  )m  +  X.V*.  (jw.  )m   .  (11) 

,  n   i  i      u  i      i   i  J  i 
i  =  0 


Since  ( j )   is  common  to  both  sides  of  (11)  and  is  not  a 


function  of  i, 


I  K 

Z   X.X*   [  Z   (-Dmd,  (jco.  )k+  (-Dkd,  (J03,  )k]  o).m 
i=0   1     k=0 


I 

=  Z   [(-l)m.  X*.V.+V*.X.]  w.m  ,  (12) 

.  n  11     ill' 

1  =  0 


where  D.  has  been  replaced  by  its  summation.   The  d,  ' s  are 
not  functions  of  i  or  m,  and 

(-1)  +(-1)   =  0,  if  m  and  k  are  not  both  even  or  odd 
=  2,  if  m  and  k  are  both  even 
=  -2  if  m  and  k  are  both  odd 


Therefore,  for  m  =  0  to  M  and  k  =  0  to  K  the  left  side  of 


(12)  can  be  expressed  as 


where 


1  =  0 


1  o  w(0,K) 

0  w(l,l)  w(l,K) 

w(2,0)        0  

0  w(3,D 

w(M,0)  w(M,l) w(M,K). 


d    (13) 


w(m,k)  =  jK  [(_l)m+(-i)K] 


d  =    d, 


(m+k) 

i 


dT 
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The  right  hand  side  of  (12)  can  be  expressed  as 


Re   (X.*V.) 

1   1 

-jIm(X1*V.)  1 
Re   (X.*V.)  .: 


(-1)M-X.*V.+V.+V.*X.   •  .M 
11111     l 


i  =  0 


(14) 


where        Re(X.*V.)    =    real   part    of   X.*V. 
11  y  11 

Im(X.*V.)    =    imaginary   part    of   X.*V. 

oo.    =    2irf.    =   

l  it 

t  =  the  time  span  over  which  v(t)  and  x(t)  are 
observed . 

Now  to  go  from  the  foregoing  general  derivation  to  the 
particular  case  where  the  order  of  the  system  is  N,  let 
M  =  K  =  N.   Then  from  (12),  (13)  and  (14)  and  after  simpli- 
fication,  the  solution  for  d  becomes 


-1  -1 
d  =  G   W  ±Z 


(15) 


where 


d  = 


do 

dl 

d2 

t 

dN 
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G 


-1 


g- 


>N 


=  1 


-T 

2i\ 


-n 


=n-2(27T') 


w 


-1 


(X.X.* 
l  l 


i  =  0 


1 

0 

w(2,0) 

0 

w(4,0) 


-1 


0      w(0,N) 

w(l,l) w(l,N) 

0      _ 

w(3,D - 

0      _ 


w(N,0)    w(N,l) w(N,N)_ 


w(m,k)  =  0,  if  m+k  is  odd 


.m+k    ,, 

i    ,  otherwise 


m,k  =  0  to  N 
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z  = 


i=0 


"zo~ 

zl 

Z2 

• 

_ZN_ 

Re(X.*V.)i   for  n  even 
1   i 


n 


-Im(X.*V. )in  for  n  odd 
1   l 

n  =  0  to  N 

The  implementation  of  the  identification  scheme  then 
takes  on  the  form  of  Figure  C2. 

For  example j  if  the  order  of  the  system  N  equals  2  and 
v(t)  and  x(t)  are  observed  through  a  time  window  of  t 
seconds  then, 


d, 


1 


0    0 

■&  ° 


/X.X.* 
/  1  1 


1  =  0 


0.   i 
.2   „ 


i2   0 


1 


2" 
4 

-  -1 

I 
/ 

/ 

i  =  0L 

Re(X.  V. ) 

l   l 

-Im(X.V  )i 
l  l 

Re(X,  *V.  )i2 

_    i   i_i 
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x(t) 
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Figure  C2.   Implementation  of 
the  Identification  Scheme 
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APPENDIX  D 
Computer  Program  PLEST  (PLANT  ESTIMATION) 

The  computer  program  is  composed  of  the  main  program 
and  subroutines  GAUSS,  DHARM,  PLANT,  RECIP  and  EA.   GAUSS 
and  DHARM  are  called  from  the  IBM/360  Scientific  Subroutine 
package  ( 360A-CM-03X)  and  are  not  included  in  the  program 
listing. 

Subroutine  GAUSS  is  used  to  obtain  a  sequence  of  real 
random  numbers  with  a  normal  distribution  and  a  selectable 
mean  and  standard  deviation.   The  expected  value  of  the 
power  spectrum  of  the  output  of  GAUSS  has  been  investigated 
by  utilizing  the  Fast  Fourier  transform  and  has  been  found 
to  approach  a  constant  value  over  all  frequencies.   GAUSS 
is  therefore  used  to  approximate  a  source  of  "white"  noise. 

Subroutine  DHARM  computes  the  Discrete  Complex  Fourier 
Transform  of  a  sequence  of  numbers.   DHARM  is  used  to 
obtain  the  frequency  spectrum  of  the  input  and  output  of 
the  system  (plant)  for  which  the  transfer  function  coef- 
ficients are  to  be  determined. 

Subroutine  RECIP  is  used  to  obtain  the  inverse  of  a 
square  matrix  with  dimensions  of  from  (1,1)  to(12,12). 

Subroutine  EA  computes  either, 


00       An 

n=0 
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or 


n 
Jo  *™T 


where   A  is  a  square  matrix 

When  the  sum  of  the  magnitudes  of  all  of  the  elements 

th 
in  the  N    term  of  (1)  or  (2)  is  less  than  or  equal  to 

-20 

10    ,  the  computation  is  stopped. 

Subroutine  PLANT  computes  the  0  and  r  matrices  of 
equation  (12)  in  Appendix  B.   The  inputs  to  PLANT  consist 
of  either  the  coefficients  of  a  transfer  function  for  a 
system  or  the  A  matrix  and  B  vector  from  a  vector  differ- 
ential equation  such  as; 

x(t)  =  Ax(t)+Bv(t)  (3) 

The  MAIN  program  (PLEST)  performs  the  normal  bookkeeping 
functions  required  of  any  comprehensive  program.   It  also 
utilizes  subroutine  PLANT  so  that  a  continuous  system  such 
as  Figure  Dl  is  converted  to  the  discrete  system  in 
Figure  D2 .   The  MAIN  program,  subroutines  DHARM  and  RECIP 
perform  the  operations  required  by  Appendix  C  to  identify 
the  system  under  test . 

The  primary  program  parameters  are  identified  in  the 
table  of  Computer  Program  Data  Cards  (Table  Dl). 
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x(t) 
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TRANSFORM 


TO  FOURIER 
TRANSFORM 


Figure  Dl.   Continuous  Version  of 
the  System  Under  Test 
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from  Equation    (12) 
Appendix   B 
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TRANSFORM 


T-N-M-K 
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Figure  D2 .   Discrete  Version  of  the  System  under  Test 
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TABLE  Dl 


COMPUTER  PROGRAM  DATA  CARDS 


Card  Format 


1    14 


2    10A8 


Name 


5    14 


14 


7   14 


NJOBS 


3   D10.4   YARN 


4   D10.4   T 


N 


M 


K 


Comments 

The  number  of  sets  of  data  cards 
to  be  operated  on. 

Arbitrary  comments  to  be  read  and 
written  on  the  print-out. 

The  variance  of  the  random  signal 
generated  by  the  subroutine  gauss 
(the  mean  has  been  preset  to  0). 

The  time  (T)  between  samples  for 
the  filter  input. 

The  number  of  successive  inputs 
to  the  filter  for  one  output. 
The  iteration  time  of  the  filter 
(TN)  is  determined  by  TN=T.N. 
0  <  N  <  12 

The  number  of  successive  inputs 
to  the  plant  under  test  for  one 
output.   The  iteration  time  of 
the  plant  is  determined  by 
TMN  =  T-M-N.   0  <  M  <  12 

th 
For  every  K    iteration  of  the 

plant  under  test  the  plant  input 

and  output  are  sampled  once  for 

the  inputs  to  the  FFT .   The  inputs 

to  the  FFT  occur  at  intervals  of 

TKMN   =    T-K-M-N.       0<K<12 
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Card  Format 
8    14 


Name 


NRUNS 


Comments 
The  number  of  times  the  plant 
coefficients  will  be  deter- 
mined.  The  resulting  coeffi- 
cients are  the  average  of 
the  results  from  NRUNS  esti- 
mates of  the  coefficients. 


14 


LTRMS 


A  lower-limit  guess  on  the 
number  of  terms  in  the  plant 
transfer  function. 
1  <  LTRMS  <  NTRMS 


10 


14 


NTRMS 


An  upper-limit  guess  on  the 
number  of  terms  in  the  plant 
transfer  function. 
1  <  NTRMS  <  12 


11 


14 


12 


14 


ISTST 


M(l) 


The  number  of  iterations  the 
filter-plant  combination  will 
go  through  to  allow  it  to 
reach  a  steady  state  condition 

2  samples  of  the  plant 
input  and  output  are  to  be 
used  on  each  run  (each  NRUNS 
from  card  8)  to  estimate  the 
transfer  function  coefficients. 

3  <  M(l)  <  10 


13 


12 


IFLT 


The  method  by  which  the  fil- 
ter will  be  characterized.  If 
IFLT  =  1  the  transfer  function 
coefficients  are  to  be  entered. 
If  IFLT  =2  the  bottom  row  of 
the  A  matrix  and  each  element 
of  the  B  vector  Is  to  be 
entered . 
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Card  Format 


Name 


14    12 


NDPFLT 


Comments 

The  exponent  of  highest  power 
of  S  in  the  denominator  of 
the  filter  transfer  function. 
0  <  NDPFLT  <  11 


15    12 


NPFLT 


The  exponent  of  highest  power 
of  S  in  the  numerator  of  the 
filter  transfer  function. 
0  <  NPFLT  <  NDPFLT 


NOTE:  if  IFLT  =  2,  go  to  (2)  16A 


(1)  16A   D10.4   DE(1) 


(1)  16B   D10.4   DE(2) 


The  coefficient  of  the  highest 
power  of  S  in  the  denominator. 


CD    16-      D10.4      DE(NDPFLT+1)      The    coefficient    of   S°    in   the 

transfer  function  denominator 


CD  17A   DlO.lJ   UN(1) 


The  coefficient  of  the  highest 
power  of  S  in  the  transfer 
function  numerator. 


(1)  17B   D10.4   UN(2) 


CD  17-   D10.4   UNCNPFLT+1)    The  coefficient  of  S°  in  the 

transfer  function  numerator. 

NOTE:   if  IFLT  =  1,  go  to  18 

(2)  16A   D10.4   A(NDPFLT,1)    The  element  A(NDPFLT,1)  in 

the  filter  A  matrix. 
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Card  Format 


Name 


Comments 


(2)  16B   D10.il   A(NDPPLT,2) 


(2)  16-   D10.4   A(NDPFLT,NDPFLT)   The  element  A  ( NDPFLT,NDPFLT) 

in  the  filter  A  matrix. 


(2)  17A   D10.4   B(l) 


The  element  B(l)  of  the 
filter  B  vector 


(2)  17B   D10.it   B(2) 


(2)  17-   D10.it   B(NDPFLT) 


The  element  B(NDPFLT)  of 
the  filter  B  vector. 


18    12      IPLT  The  method  by  which  the 

plant  will  be  character- 
ized.  Same  comments  as 
for  card  13. 

Cards  19  through  22  (for  the  plant)  are  similar  to 
cards  lit  through  17  (for  the  filter), 

NOTE:   to  omit  the  filter  enter  the  following  data 


Card 

Data 

5 

0001 

13 

01 

14 

00 

15 

00 

(1)  16A    1.0000D+00 
(1)  17A    1.0000D+00 

Repeat  Data  Cards  2  through  22  NJOBS  (from  Card  1)  times, 
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